function [final_class,final_S]=FNC_Kmeans_algorithm_PAPER(X,K,cst,dist_meas,itr_num)
% FNC_Kmeans_algorithm_PAPER groups data set X into K clusters
% using different distance measures.
%
% X: Data set to be clustered.
%    Number of rows of X = Number of genes
%    Number of columns of X = Number of experiments (cells)
% K: Number of clusters
% cst: 'a' = chi-square based distance measure (TrnasChisq, PCAChisq, PoissonC)
%      'b' = non chi-square based distance measure (Euclidean, Pearson-correlation)
% dist_meas: 'nc' = TrnasChisq
%            'pq' = PCAChisq
%            'pc' = PoissonC
%            'eu' = Euclidean
%            'neu' = Euclidean on normalized data
%            'crr' = Pearson-correlation
%            'srcrr'= Spearman Rank Correation
% itr_num: Number of iterations
%
% final_class: Cluster index for each gene
% final_S: Sum of chi-square statistic value

[num_obs,num_var]=size(X);
%% for chi-square based distance measure: TrnasChisq, PCAChisq, PoissonC
if strcmp(cst,'a')
    if strcmp(dist_meas,'nc')
        Pone=-1/sqrt(2);
        Ptwo=1/sqrt(2);
        Prest=zeros(num_var-2,1);
        k=1;
        for i=1:num_var-1
            for h=i+1:num_var
                P(:,k) = [Prest(1:i-1,:); Pone; Prest(i:h-2,:); Ptwo; Prest(h-1:end,:)];
                k=k+1;
            end
        end
        prj_X=X*P;
    elseif strcmp(dist_meas,'pq')
        covX=(X'*X)/(size(X,1)-1);
        [uu,ss,P]=svd(covX);
        prj_X=X*P;
    end

    sum_wd_old=1e10;
    for itr=1:itr_num
        % calculate 'lamda'
        randr=randperm(num_obs);
        unit=floor(num_obs/K);
        for i=1:K
            uunit=randr(unit*(i-1)+1:unit*i);
            lamda(i,:)=sum(X(uunit,:),1)./sum(sum(X(uunit,:)));
        end

        Err=1e10;
        tt=1;
        while (Err > 1e-30)
            class=[];
            for i=1:num_obs
                if strcmp(dist_meas,'pc')
                    tempX=repmat(X(i,:),K,1);
                    expt=lamda*sum(X(i,:));
                    diff=((tempX-expt).^2)./abs(expt);
                elseif strcmp(dist_meas,'nc') | strcmp(dist_meas,'pq')
                    tempX=repmat(prj_X(i,:),K,1);
                    exptt=lamda.*sum(X(i,:));
                    expt=exptt*P;
                    varc=exptt*(P.^2);
                    diff=((tempX-expt).^2)./varc;
                end
                S=sum(diff,2);
                [value(i,1),idx]=min(S);
                class(i,1)=idx;
            end

            temp_value=value;
            for h=1:K
                class_idx=[];
                temp_class=[];
                class_idx=find(abs(h-class)<1e-10);
                if length(class_idx)==0
                    [max_value,m_idx]=max(temp_value);
                    temp_class=X(m_idx,:);
                    wd(h)=inf;
                    temp_value(m_idx,:)=mean(temp_value,1);
                else
                    temp_class=X(class_idx,:);
                    wd(h)=sum(value(class_idx,1));
                end
                temp_lamda(h,:)=sum(temp_class,1)./sum(sum(temp_class));
            end
            Err=norm(lamda-temp_lamda);
            sum_wd_new=sum(wd);
            lamda=temp_lamda;
            tt=tt+1;
            if tt>5e2
                break
            end
        end
        % decide final class with minimun S
        if (Err < 1e-30) & (sum_wd_new < sum_wd_old)
            sum_wd_old=sum_wd_new;itr;
            final_class=class;
        end
    end
    final_S=sum_wd_old;

    %% for non chi-square based distance measure: euclidean & correlation
elseif strcmp(cst,'b')
    if strcmp(dist_meas,'eu')
        sum_wd_old=1e10;
    elseif strcmp(dist_meas,'neu')
        sum_wd_old=1e10;
        X=X./repmat(sum(X,2),1,size(X,2));
    elseif strcmp(dist_meas,'crr') | strcmp(dist_meas,'srcrr')
        sum_wd_old=-1e10;
    end

    for itr=1:itr_num
        centroid=X(randsample(num_obs,K),:);
        Err=1e10;
        tt=1;
        while  (Err>1e-30)
            if strcmp(dist_meas,'eu') | strcmp(dist_meas,'neu')
                for i=1:num_obs
                    dist=(repmat(X(i,:),K,1)-centroid).^2;
                    [value(i,1),idx]=min(sum(dist,2));
                    class(i,1)=idx;
                end
            elseif strcmp(dist_meas,'crr')
                Xscaled = X - repmat(mean(X,2),1,num_var);
                Xnorm = sqrt(sum(Xscaled.^2, 2));
                Cscaled = centroid - repmat(mean(centroid,2),1,num_var);
                Cnorm = sqrt(sum(Cscaled.^2, 2));
                for i=1:num_obs
                    iner_prd=Cscaled*Xscaled(i,:)';
                    mag=Cnorm.*Xnorm(i,:);
                    dist=iner_prd./mag;
                    [value(i,1),idx]=max(dist);
                    class(i,1)=idx;
                end
            elseif strcmp(dist_meas,'srcrr')
                for i=1:num_obs
                    for h=1:K
                        [dist(h,1)]=FNC_CLARITY_SpearmanRANKcorrelation(X(i,:),centroid(h,:));
                    end
                    [value(i,1),idx]=max(dist);
                    class(i,1)=idx;
                end
            end

            temp_value=value;
            for h=1:K
                class_idx=[];
                class_idx=find(abs(h-class)<1e-10);
                if length(class_idx)< 1e-10
                    if strcmp(dist_meas,'eu') | strcmp(dist_meas,'neu')
                        [max_value,m_idx]=max(temp_value);
                        wd(h)=inf;
                    elseif strcmp(dist_meas,'crr') | strcmp(dist_meas,'srcrr')
                        [max_value,m_idx]=min(temp_value);
                        wd(h)=-inf;
                    end
                    centroid_temp(h,:)=X(m_idx,:);
                    temp_value(m_idx,:)=mean(temp_value,1);
                else
                    wd(h)=sum(value(class_idx,1));
                    centroid_temp(h,:)=mean(X(class_idx,:),1);
                end
            end
            sum_wd_new=sum(wd);
            Err=norm(centroid-centroid_temp);
            centroid=centroid_temp;

            tt=tt+1;
            if tt>5e2
                break
            end
        end
        % decide final class with minimun S
        if Err<(1e-30)
            if strcmp(dist_meas,'eu') & (sum_wd_new < sum_wd_old)
                sum_wd_old=sum_wd_new;itr;
                final_class=class;
            elseif strcmp(dist_meas,'neu') & (sum_wd_new < sum_wd_old)
                sum_wd_old=sum_wd_new;itr;
                final_class=class;
            elseif strcmp(dist_meas,'crr') & (sum_wd_new > sum_wd_old)
                sum_wd_old=sum_wd_new;itr;
                final_class=class;
            elseif strcmp(dist_meas,'srcrr') & (sum_wd_new > sum_wd_old)
                sum_wd_old=sum_wd_new;
                final_class=class;
            end
        end
    end
end
final_S=sum_wd_old;
end